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A Multifunction Space-Borne Scanner 
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in Chinese No 2, 30 Apr 90 pp 52-56 


[Article by Kong Huixing [7895 1920 5281], Shanghai Institute of Technical 
Physics] “ 


[Text] Abstract 


A proposal is made to add five detection bands to an existing very-high- 
resolution (VHR) scanning radiometer to enhance its capability for the 
observation of meteorological phenomena, ocean environment, and natural 
disasters. 


(1) Performance of the Existing Very-High-Resolution Scannin, Radiometer 


Test results show that the performance indices of China's existing VHR 
scanning radiometer in the visible-light and infrared band meet all the 
design requirements, as shown in Table 1. Both visible-light and infrared 
channels produce clear images, and have good observation capability for 
cloud maps, ground covers, lakes and ocean water. 


(Il) Satellite and Radiometer Design Requirements for Monitoring the 
Environment and Natural Disasters 


Because natural disasters occur abruptly in time and in space, real-time 
monitoring is required in order to provide timely reporting and assistance. 
An effective approach of monitoring the entire Chinese landmass is to use a 
satellite-borne VHR scanning radiometer operating at high repetition 
frequency. 


The results of recent studies in earth science and environmental science 
show that in order to explain the evolution of certain important environ- 
mental phenomena, it is necessary to consider the earth as a single entity, 
and to study the global effects between the atmosphere, the ocean, and the 
land, the gain and loss of earth radiation, and the combustion of fossil 
fuels, etc. To accomplish this requires the accumulation of over 10 years 
of observations with high spatial and time resolution, as indicated in 
Table 2. 








Table 1. On-Orbit Performance of Very-High-Resolution Scanning 





Radiometer 








Parameter Design specification Test results 
Detection 0.58-0.68 um Visible cloud maps, ground cover, cloud 
band and snow separation, ocean sediment 
0.725-1.1 um Visible cloud maps, ground cover, cloud 
and snow separation, water and land 
border separation 
0.48-0.53 um Ocean water color, chlorophyll 
| distribution 
0.53-0.58 um Ocean water color, chlorophyll 
distribution 
10.5-12.5 um Infrared cloud maps, water and land 
border separation, temperature 
measurement 
Spatial 1.08 km 1.08 km 
resolution 
Detection Ch,, Chg, NEAp<0.16% NEAp<0.11% 
sensitivity Ch3, Chg, NEAp<0.25% NEAp <0.15% 
Chs NEAT<0.8 K NEAT<0.5 K 
Detection Ch3,Chy, +10% | +10% 
accuracy Che 1.0 K approximately 1.0 K 





Table 2. Time and Spatial Resolution Requirements for Environmental 
Satellite Observations 








Object of | Ground 
observation | Repetition period resolution 
Atmosphere At least 0.5 day/observation 1 km 
Ocean 0.5-1 day/observation (for ocean conditions) 0.5-1 km 

3-5 days/observation (for medium and large-scale 
| ocean phenomena) 1-4 km 
Ground cover | 7-10 days/observation 0.5-1 km 








Although earth resource satellites have very high spatial resolution, it 


cannot satisfy the repetition-frequency requirement for environmental 
monitoring because it takes 16-26 days to cover the whole earth. 
tion, the multi-band scanner of an environmental satellite has 5-10 times 
higher sensitivity than that of an earth resource satellite (see Table 3) and 


therefore can detect minute changes in the environmental parameters and 


predict its future development. Because of these unique features of the 


environmental satellite, it should be further developed to provide the 
capability for monitoring the environment and natural disasters. 
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In addi- 
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(III) 10-Band Scanning Radiometer 
In order to enhance the monitoring capability of the existing VHR scanning 
radiometer, it is necessary to add five more detection bands to the existing 
bands to provide a 10-band observation capability. Specifically, the 
enhancements include the following: 


(1) Adding a 3.55-3.93 um medium-infrared band for monitoring forest fires 
and grassland fires; the low water absorption in this band will provide 
improved temperature measurement accuracy over the ocean surface in the mid- 
low latitude region. 


This band has very high sensitivity with respect to fire because it is 
located at the radiation peak of a 780 K target. The following table shows 
the calculated results of radiation luminance for targets at different 
temperatures. 


Table 4. Radiation Luminance of Targets at Different Temperatures in the 
10.5-12.5 um Band and the 3.55-3.93 um Band 





Target temperature | 
~ T=1100 K | T=800 K = T=300 K__N(1100 K) | N(800 K) 


Target luainance ~~~ N(300 K) , N(300 K) 





abe + —— 
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3.55-3.93 um | 
Luminance (W/cm2*ster) 


10.5-12.5 um 
Luminance (W/cm2*ster) 


0.1924 5.08x1072 -1.7x107>2 1.13x104  3.9x103 


1 5.59x10-2 |3.13x1072 1.86x1073 30.0 | 16.8 
| 





The radiated power received by a scanner from a target which does not fill 
the instantaneous field of view is proportional to the target area. Thus, 
the radiated power of a 20x20 m2 fire at 800°K is equivalent to that of a 
1.1xl.1 Km target at 300°K. Taking into account atmospheric attenuation 
and margin threshold setting, it is estimated that the scanner can detect a 
30x30 m2 open fire and issue a fire alarm. 


(2) Adding a 1.58-1.64 um shortwave infrared band for measuring the water 
content in soil and crops and for monitoring drought conditions. 


The soil reflectivity drops sharply in the 1.3-2.5 um band when it contains 
water; in particular, there are water absorption peaks at 1.44 um, 1.90 um 
and 2.2 um. When the soil moisture is in the 0.8 percent to 20.2 percent 
range, an increase of 1 percent in water content will reduce the reflec- 
tivity by 1 percent (see Figure 1). The water content in the leaves of 
ground vegetation will clearly affect the reflectivity in the 1.3-2.5 um 
band (see Figure 2). The sensitivity of the scanning radiometer is estimated 
to NEAp<0.1%, which is sufficient to detect a 1 percent change in water 
content; therefore, the scanner is effective in monitoring drought condi- 
tions. Considering the fact that there are H70 and CO) absorption in the 
atmosphere, the detection band is chosen to be 1.58-1.64 um. 
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Figure 1. Spectral Characteristics of Soil Reflectivity With 
Different Moisture Content 
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Figure 2. Reflectivity Spectrum of Clay and Wheat Leaves With 
Different Moisture Content 


(3) Adding a 0.433-0.453 wm band for measuring the chlorophyll content in 
the ocean. This band supplements the existing 0.48-0.53 um band and the 
0.53-0.58 wm band to measure the water color which reflects a wide range of 
chlorophyll content. The collected data can not only be used for preliminary 
survey of ocean productivity but also be used in conjunction with surface 
temperature data for fishing-ground forecast. 


(4) Because of scattering by molecules anc colloidal particles in the 
atmosphere in the 0.433-0.453 um band, the sky background is one order of 
magnitude brighter than the ocean surface. In order to isolate the water- 
color signal of the ocean from the signal of the bright sky background, it 
is necessary to add a 0.745-0.785 um band to correct for the effect of che 
colloidal particles in the atmosphere. 


(5) The existing 10.5-12.5 um thermal infrared band of the scanning 
radiometer is split into two bands: 10.5-11.5 um and 11.5-12.5 um, in order 
to partially cancel the atmospheric effect so that higher measurement 
accuracy of the ocean temperature (approximately 1.5 K) can be achieved. 


(6) Because of the large differences between the reflectivities of soil, 
ground cover, and water in the 0.725-1.1 um band, it is an important band 
for monitoring land-water separation, river flow rate, and flood damage. 








In order to improve the capability of separating water and land, the ground 
resolution in this band is increased to 540 m. 


(7) 


Currently, the 0.725-1.1 wm band and the 0.58-0.68 um band are combined 


to monitor the ground cover index, which is a measure of ground cover 
distribution, growth trend of crops, and estimated wheat production. In 
order to match the spatial resolution between these two bands, the resolution 


in the 0.58-0.68 um band is also increased to 540 m. 
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Figure 3. Upward Spectral Luminance of Ocean Water With 


Different Chlorophyll Content 
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Figure 4. Reflectivity Spectrum of Typical Targets 


The modified scanning radiometer has 10 detection bands. It includes the 
existing five bands of the AVHRR of the NOAA weather satellite and the 
1.58-1.64 um band which is to be added on the NOAA-K, L, M satellites. 

The addition of the 0.433-0.453 um band and the 0.745-0.785 um band to the 
existing two ocean channels of the scanning radiometer provides the quanti- 
tative measurement capability of chlorophyll content and sediment level in 
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the ocean. 


Therefore, the improved 10-band scanning radiometer will have an 


integral capability of monitoring meteorological phenomena, water color as 
well as natural disasters such as floods, fires, and drought. 


The author wishes to express his thanks to Mr Kuang Dingpo and Mr Zheng Qinpo 


for their assistance in this work. 


Table 5. Performance of the 10-Band Scanning Radiometer 














Channel | | Spatial | Coverage | 
No Band | resolution | period Application 
| 
1 0.433-0.453 um | 1.08x1.08 km2 | 24 hours | Ocean chlorophyll 
2 0.48-0.53 um | 1.08x1.08 km? | 24 hours | Ocean chlorophyll 
3 0.53-0.58 um 1.08x1.08 km2 | 24 hours Ocean chlorophyll, yellow 
| | | substance 
4 0.58-0.68 um = 0.54x0.54 km? | 24 hours | Cloud maps, ground cover, 
| ocean sediment 
5 -0.725-1.1 um | 0.54x0.54 km2 24 hours Cloud maps, ground cover, 
water and land separation 
6 0.745-0.785 um | 1.08x1.08 km2 , 24 hours | Atmospheric correction 
7 | 1.58-1.64 um /1.08x1.08 km2 , 24 hours | Drought detection, cloud 
| | | and snow separation 
8 3.55-3.93 um | 1.08x1.08 km2 12 hours | Fire detection, ocean 
! | surface temperature 
9 10.5-11.5 um |1.08x1.08 km2 , 12 hours | Infrared cloud maps, 
| | | temperature measurement 
10 ; 11.5-12.5 um | 1.08x1.08 km2 12 hours | Infrared cloud maps, 
| | | | temperature measurement 
References 

1. Zheng Qinpo, et al., “Weather Satellite Remote Sensor and Evaluation 
of the Quality of Cloud Maps," 1988, 10. 

2. Fang Zhongyi, “Application of Weather Satellite Data," 1989, 8. 

3. Gordon, H. R., Clark, D. K., Muller, J. L., and Hovis, W. A., 
"Phytoplankton Pigments From the Nioubus-7 Coastal Zone Color Scanner: 
Comparisons With Surface Measurements,'' SCIENCE, 210, 63, 1980. 

4, Esaias, W., et al., "Moderate Resolution Imaging Spectrometer," NASA EOS. 
1986. 

5. "Multispectral Scanning System and Their Potential Application to Earth 


Resource,’ SURVEY, 1972. 








AEROSPACE 


Dynamics of Simplified Model of Solar Panel and Supporting Components 


90FEO179B Beijing YUHANG XUEBAO [JOURNAL OF CHINESE SOCIETY OF ASTRONAUTICS] 
in Chinese No 2, 30 Apr 90 pp 57-64 


[Article by Zhang Jiazhong [1728 0857 6988] and Shao Chengxun [6730 2052 
8133], Harbin Polytechnical University] 


[Text] Abstract 


Expcrimental results and theoretical analysis of a simplified model of a 
solar panel and its supporting components used on a space station are 
presented. Specifically, the dynamic characteristics of the system when it 
is subject to disturbances such as attitude adjustment and orbit maneuvers 
are investigated. The results show that because of the inertial non-linear 
coupling effect of a large flexible space structure, the presence of 
internal resonance will produce abrupt changes in the dynamic response of 
the structure when the amplitude of the disturbance torque exceeds a 
threshold value. These changes exhibit non-linear characteristics of 
simultaneous excitation of the coupled modes. The abrupt increases in the 
amplitudes are expected to cause difficulties in maintaining the shape and 
controlling vibration of the flexible structure. 


I. Introduction 


For a large space structure such as the space station which operates in 
near-zero load conditions, the strength requirements are very low. To 
reduce the weight of the structure, many of its components are extremely 
flexible. However, these components are subject to large overload during 
orbit injection; therefore they are compactly stowed during launch, then 
deployed or re-assembled after achieving orbit. A typical example of such 
flexible component is the foldeble solar panel discussed in reference l. 
Figure 1 shows the configuration of a postulated space station, which con- 
sists of the extremely flexible solar arrays, the radiation panels and a 
large antenna. 


Because a space station generally undergoes a series of on-orbit operations 
such as deployment and re-assembly of its components, docking with the 
spacecraft, attitude control and orbit maneuvers, it is inevitably subject 
to different types of disturbances. A large flexible structure such as the 











space station has a high concentration of low-frequency natural modes. An 
analysis of a postulated space station? shows that there are 14 natural modes 
in the frequency range 0.21-0.75 Hz, and 20 natural modes within a 2-Hz 
range. These include the modes of the flexible structure as well as the 
aggregate modes of the space station. Furthermore, these flexible components 
which are connected through the main body of the space station are mutually 
coupled. The non-linear coupling has a significant effect on the dynamic 
characteristics of the system. Reference 3 pointed out that these non-linear 
coupling effects will cause local vibrations in the structure. References 

4 and 5 have analyzed the effect of inertial non-linear coupling on the 
dynamic characteristics of a flexible structure, and pointed out that the 
presence of internal resonance in the system will cause abrupt changes in 

the system dynamic response even if the coupling is weak. For a large 
flexible structure (Figure 1) with a large number of natural modes, the 
presence of internal resonance is almost unavoidable. Therefore, this 

effect must be carefully studied. 














Figure 1. Proposed Space Station 


The dynamic equations of a large space structure expressed in modal coordi- 
nates are of the following form: 


~ 
i 


, 2 Q i~ . — . ; ; © 2 
Gp +079 = Fy | Dada (Avigidit BeiigsGi) t2@rsrGe | 
+eBrd order terms] (1) 


where e€ is a small parameter, and the second term on the right side of the 
equation is the second-order inertial coupling term. Q,; is the generalized 
force, My is the generalized mass, wy is the frequency of the r-th order 
linear mode, and A;yj; and Byj4 are the corresponding coupling coefficients. 
Under normal conditions, the effect of these non-linear coupling terms can 
be neglected, which leads to the conventional linear dynamic equations: 


2 Q, 
G-+2a@0-G-+@-9e = M. (2) 








Equation (2) provides an accurate description of the system characteristics 
and yields the correct system response. But when internal resonance exists 
between the natural modes of the system, i.e., when the relationship 


Wi =O; (3) 
exists between two natural modes, or when the relationship 
Dk= @i TO (4) 


exists between three natural modes, or when both relationships are present, 
then as the amplitude of the disturbance increases, abrupt changes of the 
system dynamic response may occur. In this case, the effect of non-linear 
coupling can no longer be neglected. This article presents the experi- 
mental results and theoretical analysis of the phenomenon of internal 
resonance under non-linear coupling for a simplified model of solar panel 
and its supporting components. 


Il. Experiments 


In order to study the effect of inertial non-linear coupling on svstem 
dynamics, we have conducted an experimental investigation and performed 
theoretical analysis of a simple model of solar panel and its supporting 
system as shown in Figure 2. The solar panel is modeled as a beam under 
bending and torsional loads, the supporting system is modeled as another 
beam. We only consider the case where the solar panel is perpendicular to 
the supporting beam; the physical parameters of the two components are 
adjusted so that the natural bending mode in the plane of the solar panel 
and the supporting components is w); = 6.1 Hz, the natural torsional mode 
about the axis of the supporting beam is wr = 56 Hz, and the natural 
torsional mode of the solar panel about its own axis is w,. = 52 Hz. Clearly, 
these three frequencies satisfy the following internal resonance relation- 
ship: 


Wr ==, t+ O¢ (5) 








Figure 2 
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A sine-wave excitation signal is generated using a SCO200 scan-controlled 
signal generator through a LDS resonator; its spectrum and time history are 
recorded using a spectrum analyzer and a X-Y recorder. The torsional vibra- 
tion of the supporting beam is measured by a micro-accelerometer located at 
the end of the supporting beam. The torsion and bending of the solar panel 
are measured by the strain gauge located on both sides of the solar panel. 


Figure 3 shows the torsional frequency response of the supporting beam under 
a sine-wave scan. Figure (a) shows the response under low excitation level, 
which is clearly linear and has a peak response at 56 Hz. Figures (b), (c) 
show that as the excitation level increases, the non-linear characteristics 
become increasingly apparent, and multiple peaks appear on the frequency 
response curve. Figure (d) shows the typical non-linear behavior of the 
frequency response as the scan direction is reversed. 
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Figure 3 


Figure 4 shows the torsional spectrum of the supporting beam, and the bending 
spectrum and torsional spectrum of the solar panel under a 56-Hz sine-wave 
excitation. It can be seen that at low excitation levels (Figure (a)), only 
torsional vibration is produced; at high excitation levels (Figure (b)), a 
drive frequency of wy will excite not only torsional vibration of supporting 
beam, but also coupled bending and torsional vibrations of the solar panel. 
In other words, at high excitation levels, the non-linear effects become very 
pronounced; the torsional disturbances acting on the supporting beam are 
transmitted to the solar panel through internal resonance, producing bending 
and torsional vibrations in the panel. 
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Figure 5 illustrates the effect of the disturbance level: curve 1 is the 
torsional vibration of the supporting beam, and curve 2 is the vibration of 
the solar panel; the horizontal axis denotes the excitation level. These 
curves clearly show that when the excitation level is less than a certain 
critical value, the system dynamics remain in the linear domain, where 
torsional vibration is produced only in the supporting beam, the solar panel 
is free of vibration. When the excitation level is greater than this 
critical value, the system dynamics become non-linear; the disturbance energy 
acting on the supporting beam is transmitted to the solar panel through 
internal resonance, and produces vibration in the panel. As the excitation 
level further increases, the amplitude of the torsional vibration of the 
supporting system reaches a saturation level, whereas the amplitude of vibra- 
tion of the solar panel continues to increase. 
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Figure 5 
For a large flexible structure with numerous low-frequency natural modes, 


this internal resonance phenomenon clearly is a potentially dangerous 
problem. 
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III. Theoretical Analysis 


The solar panel and its supporting system are represented by the simplified 
model shown in Figure 6. Assume that the disturbance torque is about the 
axis of the supporting beam and has a sine-wave excitation with frequency 2”. 
Let the supporting system and the solar panel be represented by beam 1 and 
beam 2 respectively, both of which are subject to bending and torsion. 

















Figure 6 


The relative deformation of beam 2, denoted by u(z,t) and w(z,t), can be 
expressed as follows: 


u(z,t) =u.f,(2) (6) 
py (z+t) = gog(z) (7) 

where 
Uy > Ww, Ko! (8) 


Vol is a component of the eigenvector of the first-order bending mode of the 
system, and W; is the generalized coordinate of the bending mode in the plane. 
Yo is the coordinate of the torsional mode about the beam 2 axis. Let ao be 
the generalized coordinate of the system torsional mode about the beam 1 axis, 
then the system dynamic equations expressed in modal coordinates are: 


. ; T, VM ip as a +- 
Qo 7A = Vo J. cos Qt — 261ra: a) — Lp Bis lH ot 2h Poti Po (9) 
ae My ee, 
ot OePfo= —26becaePo a “F  PoBigaoll 1 (10) 
M 


W.+o0iW,= —26,0,W,- VY iPorBirdinGo (11) 


where Bj is the second-order coupling coefficient: 
7 
Bie=\ (1-2) fi(2) (2) dz (12) 
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It can be seen from equations (9)-(11) that the coupling term [Wi¢o + Wig + 
W199] on the right side of equation (9) absorbs part of the direct excitation 
of the ao response; the coupling terms GoW) and Gogo on the right side of 
equation (10) and equation (11) cause indirect excitation of the torsional 
deformation about the beam 2 axis #9 and of the bending deformation Wy 
respectively. These coupling terms are important if the following internal 
resonance condition exists: 


(W 7 @, Tay (13) 


Otherwise they can be neglected. This conclusion and the experimental results 
presented in the previous section can be verified by theoretical analysis of 
this simplified model using the multiple scale method. 


To proceed, we introduce an arbitrary small parameter ce and re-write 
equations (9)-(11) in the form: 


Y +07), =e[reosQt—2trarY ,— uk, (X X14 2K,N.4+4 KX, X,)) (14) 
\, + G); X | é| —2é:0@ X, ~ iN a’ XxX.) (15) 
V,taiXx, =e( — 2b eae X,-Kiu¥.X J (16) 
where 
} Qos el X nH, el X r¢ 
Af r l. K - ty B 

M, 2 2 Tl, ., ats fe) 

a ae adel Mi ie (17) 


r is the radius of gyration, M, is the generalized mass corresponding to the 
system bending mode, Ip is the polar moment of inertia of beam 2, My is the 
mass at the end of beam 2, To is the amplitude of the applied torque, Jy) is 
the generalized rotational inertia, yo is a component of the torsional mode 
at the point of the applied torque, and ? is the length of beam 2. 


By introducing a slowly time-varying parameter Tj = et, the solution of 
equations (14)-(16) can be expressed in the form: 


VY. (t) =6,(7 ,) coslowrt+P.(T,)] +eY,,(:,T7:) (18) 
Y,(t) =a,(T,) cosLo,t+a,(T,)] +eX,,(57,) (19) 
Y,(t) a,(7T,) cos(met+a,(T,)]+eX;,(t,7,) (20) 


Equations (18)-(20) represent the torsion of beam 1, the bending of the 
system in the plane, and the torsion of beam 2 respectively. It can be seen 
from these equations that each solution corresponds to a domain of its own 
natural frequency. From the bounded conditions of the disturbance term 
Y11(t, Tj), ome can derive a first-order differential equation for each of 
the response amplitudes bj, a}, and ag with respect to the slowly varying time 
Tj. The degree of resonance can be defined by two parameters o and 0o9: 
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$2 WrwTreé (21) 
Or= Or tert ea, (22) 


Expanding the differential equations for the modal responses aj, a9, and bj, 
one obtains: 


, r j q 7 | rms 

@,a 4: @;da; ’ K 6 g:-0:51N YY; (23) 
, l } . . 
@,a:a; 4 \ a2 1:0-COS/; (24) 
Wea). Crea: t oF K ya’ a,0 sin y (25) 
] 

Wea: om IK @Wra,;0;COS) (26) 
Orb ~ ,@'56 “4 uk (@ Vg ] ny — r (27) 

; l 
'@) b6,f 4 uly \¢ & Bas Cl r > f (28) 

where 

f O p (29) 
¥.=0,T,-a.-@: +6 (30; 


i. 


In the above equation, the superscript denotes derivative with respect to 
T,- The system equations (23)-(28) have two steady-state solutions. One of 
the solutions is: 


ai 0. ]: ' (31) 
b= so aEA (32) 


which indicates that the indirectly excited motions are stationary, and the 
response of the directly excited motions is linear. 


The second steady-state solution is: 


. , - © | y +- (7 ) [ 
b, = $919 ae [14 _ | (33) 
7 \ i? . ,W ~Ge@Me): J 
; E10: 
Qa: = - 2 a; (34) 
GeWe 
9 . yf 2 . ] 
— ——— - | 2. Pytv rr (35) 
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where 


- lsCE, wr (E,0, + Fede) Oo\lo +o,) ] (36) 

l,=I0(€.0,+§<@0e)0+é7a@7(0+0,)] (37) 
r. 8@:0oV Siée 

( rk 12(4¢,@ Tv EeWe) (38) 


This solution represents the non-linear response of the system. The three 
types of motions are excited simultaneously; also, when o = 0], the 

torsional amplitude b), has a minimum value which is independent of the level 
of excitation. This is the so-called saturation phenomenon which is clearly 
evident in the experimental results of Figure 5. The amplitudes of a) and a? 
increase with increasing excitation level. This phenomenon can also be 
observed in Figure 5. 


Equation (35) also shows that whether or not a stable solution exists depends 
on the system parameters. 


The to value of the critical point can be obtained from equations (32) and 
(33): 


*. \ a ae K.. i (f,0, 4 EeWe) (39) 
Clearly, to is a minimum when 0, +08 0 and o = 0: 
Come =Be,00E7_[8!°° (40) 
V . 


However, 0] + o = O and o = O imply that the internal resonance relationship 


@, TDs 
and the resonance relationship { = wy exist. 


Figure 7 shows the theoretical curves of the amplitudes of the modal 
responses as a function of the excitation parameter t. When Tt exceeds the 
critical value, the amplitudes of the coupled modal responses vary 
differently: one approaches saturation, and the other two undergo abrupt 
changes, then continue to rise sharply. 


IV. Conclusion 


In this article, the effect of inertial non-linear coupling on the dynamic 
characteristics of a flexible system is studied using a simplified model of 
solar panel and its supporting system. Experimental results and theoretical 
analysis show that when the external excitation is in resonance with one of 
the non-linear coupled modes, and internal resonance also exists between the 
coupled modes, the non-linear coupling can produce abrupt changes in the 
dynamic characteristics of the system. The effect of internal resonance 
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produces indirect excitation of the coupled modes which are independent of 
the external excitation. For a large, flexible structure such as the space 
station, there are numerous sources of excitation and a high concentration of 
natural modes which may cause internal resonance. This type of internal 
resonance is a potential threat to vibration suppression, configuration 
maintenance, pointing accuracy and control accuracy for a large flexible 
structure. Therefore, detailed investigation of this effect of non-linear 
coupled vibration is essential to the design of space station. 
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[Text] Abstract 


A method of coarse geometric rectification of full-size satellite panoramic 
images is proposed. The method is based on certain intrinsic geometric 
relationships between the satellite image and the ground, and takes into 
account the physical process of satellite imaging. The rectification of 
satellite images is accomplished by applying parameter optimization technique 
which does not depend solely on the satellite parameters, but also relies 

on a small number of carefully selected and accurately measured ground 
control points as well as a geometric model. Simulation results show that 
the proposed method is very effective for the coarse processing of satellite 
panoramic images. 


I. Introduction 


Ground observation and imaging by a retrievable remote-sensing satellite is 
accomplished by panoramic scanning. Because of the low satellite orbit and 
the large scan angle, the images suffer severe geometric distortions, 
particularly at the edges where errors can be as large as several 
kilometers. In order to promote and apply the satellite panoramic imaging 
technique, it is essential to find a practical method of geometric rectifi- 
cation. 


The geometric methods which are currently used for satellite image 
processing include the model method and the polynomial method.1 The model 
method involves establishing a mathematical model to describe the satellite 
imaging process and requires accurate knowledge of the satellite orbit, 
satellite attitude and camera parameters; the accuracy of these parameters 
will have a significant effect on the accuracy of geometric processing. The 
polynomial method does not require knowledge of satellite parameters, but it 
does require a large number of accurately measured and uniformly distributed 
ground control points; this method produces rather poor results at the edges 
of the images. 


18 











The method proposed in this paper is based on the geometric relationships of 
satellite imaging and tries to avcid the deficiencies of the above two 
methods. Specifically, it takes into account the physical process of 
satellite imaging, but does not rely solely on the satellite parameters; 
rather, it uses the intrinsic geometric relationship between the satellite 
image and a small number of accurately measured ground control points and 
performs geometric rectification of the satellite image by applying optimal 
parameter estimation technique. 


II. Mathematical Model2 


To simplify the analysis, the earth is assumed to be a sphere of radius R, 
and the altitude of the satellite is denoted by H. Since the effective area 
of each film is approximately 200x1400 mm2, its width corresponds to a size 
of only several tens of kilometers on the ground. Hence, the projected scan 
lines can be approximated by straight lines, and the ground segment corres- 
ponding to each film can be approximated by a curved surface taken from a 
cylinder of radius R, as indicated by the surface EFPB in Figure 1. Clearly, 
for a scan angle 8, one can apply the law of sine to the triangle AOSP, to 
obtain the following equation: 





2a= sin~'[ sin p - ( - ~) I-A (1) 
R- sin2 
ary a (2) 


Recognizing that AOS,P, is an isosceles triangle, one can determine the 
corresponding angles, and apply the law of sine to the triangle ASS)P, to 
obtain: 


|= H cosa/cos (a+ Bf) (3) 





Figure 1. Geometric Relationship of Satellite Imaging 
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1. Geometric Relationship Along the Scan Direction 


As shown in Figure 2, in the ground coordinate system 0-XY, Lo, is the ground 
projection of the panoramic scan slit when the scan angle 8 = 0; it is repre- 
sented by the equation Y = KX+E. At time tj, the ground projection of the 
slit is Lj, and an image of the straight line L;' which passes through the 
point (Xj, Yj) is taken. The distance between Lj and Lo is vtj, where v is 
the translational velocity of L, due to earth rotation, which is to be deter- 
mined. It is a function of the inclination of the satellite orbit, the 
satellite velocity, the earth rotational speed and the latitude and longitude 
of the satellite subpoint. Dy is the distance between L;' and Lj. 
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Figure 2. Geometric Relationship in the Ground Coordinate System 


In the original film coordinate system O-xy shown ii. Figure 3, it is assumed 
that the satellite attitude remains unchanged during the imaging process; in 
other words, during this period the equation of the theoretical center edge- 
contact line (8 = 0) remains constant 25: y = kxte. By "theoretical" it is 
meant that % passes through the primary image point, and the image line at 
any time is parallel to 2). dj is the distance between 2; and £,, and 
corresponds to the ground distance Dj. 








Figure 3. Geometric Relationship in the Film Coordinate System 


From the above analysis, one can readily derive: 


_IKXi+E-Y;i| = (4) 


ad J/1+K°? 
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lkxi; te—yi| 
—> Gowers , 5 
“ J 1+k’ (>) 


where K, E, k, e, v are to be determined, ty = dy/(v, ¥ 1+k2), v, is the film 
speed of the camera. 


Substituting 8; = d;/f (f is the focal length of the camera) into equation 
(1), one can derive the following relationship between D, and dj;: 


cr \ 
D;=R -{ sin-f sin 9414 fale a (6) 
It is clear that even if errors are present in the satellite actitude and 
orbit parameters, as long as they remain unchanged during the imaging 
process, this geometric relationship also remains unchanged. In other words, 
equation (6) shows that the accuracy of geometric processing is independent ~ 
of errors in the satellite attitude parameters, which is a major advantage 
of this method. Since each film only takes 3 seconds to complete during the 
imaging process, it is reasonable to assume that the effective parameters 
remain unchanged during this period. 


2. Geometric Relationship Perpendicular to the Scan Direction 


To account for small variations in satellite motion, earth rotation and the 
direction of film advancement, we use a second-order curve to represent the 
scan line of the slit-center and its ground track; these equations are of 
the form: x = ay2+by+c (film) and Y = AX2+BX+C (ground), as indicated in 
Figure 2 and Figure 3. In the ground coordinate system, let us draw a 
straight line from the ground control point (Xj, Yj): (Y-Yy) = K(X-Xj), 

which intersects the center scan line Y = AX2+BX+C at (X,, Y,.), then we have: 


DY = S(X;—X)?*+(Vi- Ye? (7) 
Similarly, in the original film coordinate system, we draw a straight line 


from (xj, yi): y-yy = k(x-xj), which intersects the center scan line 
x = ay2+bytc at (x5, yo), then we have: 





dy =v (x; =x)? = (¥i— Ys)? (8) 


Now, we shall derive a geometric relationship between DYy and dy;. Figure 4 
shows the geometry between DY; and dy; where a pitch angle ¢ is present. 

8 is the angle of the scan beam. It should be pointed out that since y and 
dy can be either positive or negative, they can be combined to form several 
different geometries; however, the derivation in this paper is based only on 
the geometry shown in Figure 4; it can be readily extended to other 
geometries. From Figure 4, it is clear that 


el 








DY =! - (tg(@+6) —lgq] 


sin2¢" tgg + lyd 


"K sing .1- temtgd ~ eg 7 


cosa r ftg¢g+dy 7 


7 to a SP | 
cos a+ 4 )L f—teedy 18" J (9) 


a fy . 


Equation (9) is affected only by the pitch angle ¥; it is independent of other 
attitude parameters. 








== _ 
Figure 4. Geometric Relationship Between DY and dy 
3. Objective Function 
In this paper, the physical parameters are cetermined by solving the 


parameter optimization problem using the Powell method. The objective func- 
tions to be optimized are: 


Pi (hyve, AVE Vv) >"! D, Ke sin =| sin < Hite) |G / (10) 
filavbsc, AVBV Cy) S DY.-H — 
| eoslat %) 
it di, 2 
t= ~ -dy, 7 leg (11) 


4. Geometric Rectification 


Having determined the parameters {K, E, k, e, v, A, B, C, a, b, c, #} using 
the Powell method, the process can be reversed to establish the relation- 
ships between corresponding points on the original film and the new film 
after rectification. A difficult problem with this process is establishing 
the quadrant relationship between corresponding points in the new picture and 
the original picture. 
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Consider a point p(x,y) on the original film and its corresponding point on 
the new film P(X,Y). The computational procedure and the key formulas are 
given below: 


1) Determine xo; Yo 


x,=av*,-9ON 
2) Calculate d, dy 
by an 
d 
v 1+R 
dv _ y / Y —x, k | y as v)- 


3) Calculate D, DY from equations (6) and (9) 


4) Solve the simultaneous equations 


Yy -Y, ~ KR (¥ -X,) 

Y= AX 424,40 

DY =J/(X -%,)°=(V—¥,) eX Fy XY) 
D- AY+F—-] - ee 


‘\ l s a5 


This completes the transformation between the points p(x,y) and P(X,Y). It 
should be noted that the positional relationship between corresponding points 
on the original film and the new film is controlled by the four "quadrants" 
formed by the center edge contact and the center scan line; i.e., the 
"quadrant" relationship between corresponding points remains unchanged. 


To save computer time in the actual rectification process, we only use the 
above model to calculate the coarse lattice points; then we apply the finite 
element method (first-order approximation inside a triangle) to calculate 
points within a lattice. 


5. Other Considerations 


Since the panoramic scan camera uses mechanical drive mechanisms for film 
advancement, jitters will inevitably appear, causing errors in the images. 

To solve this problem, we introduce in the objective function an additional 
parameter 6, which is a function of d;. The form of this function is shown 
in Figure 5, where A is the residue at d;, and L is the range of influence 

of jitter, which is the parameter to be optimized. For a control point whose 
residue is larger than a certain value, optimization is performed on an 
objective function with the added parameter 6(dj), which determines a model 
function that reflects the jitter of the entire frame. 
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Figure 5. Simulated Jitter Function of Film Advancement 


In addition, because of variations in satellite attitude parameters during 
the scan process and because of dimensional considerations, the film 
digitization process must be carried out frame by frame (due to the limita- 
tion of C-4500). Therefore, the center scan line or its ground track cannot 
be modeied by a single second-order curve; it must be modeled within each 
frame to achieve the required processing accuracy. 


Another factor that should be considered is the effect of difference in 
surface altitude. This effect is particularly pronounced for domestic 
satellites as they are in low-earth orbits and have large scan angles. In 
this article, the effect of difference in surface altitude is taken into 
account by introducing a correction factor to the orbit altitude H in the 
geometric relationship. 


III. Test Results 


In order to demonstrate the feasibility of the proposed method, we have 
conducted test on a full-size satellite panoramic image. The dimensions of 
the image are approximately 200x1400 mm2. The image is divided into six 
frames, each of which is digitized using the C-4500; the individual frames 
are joined together in the computer using identical points to form a com- 
plete digitized image. The control points are selected by performing a 
correlation analysis between the image and the terrain map. However, 
because a large portion of the area covered by this image is desert, 
selection of control points is very difficult, and they are not uniformly 
distributed. As a result, the RMS processing errors using the polynomial 
method or the traditional model method are of the order of 800-1000 n, 

which clearly do not meet the accuracy requirements. In this paper, a total 
of 50 ground target points are selected, 30 of which are control points, and 
the remaining 20 are check points. After processing the data using the 
method described in this paper, the RMS error is o = 206 m, which fully 
meets the accuracy requirement of coarse geometric processing for satellite 
images. 


This method can also be used for precision rectification of sub-images of a 
region with a large number of uniformly distributed control points. For 
example, a test was performed on a 110x110 mm2 sub-region using 30 control 
points; the RMS error of the rectified image is o = 45 m. 
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Higher processing accuracy can be achieved for a sub-image because: 1) there 
are larger number of control points; 2) the various parameters remain essen- 
tially constant within a shorter time interval; and 3) it is easier to model 
the center scan line and the center edge contact. However, precision 
processing the full-scale image all at once would be difficult. 


IV. Conclusion 


The avove test results show that the proposed method is effective for coarse 
geometric processing of satellite panoramic images. This mthod is particu- 
larly useful in cases where the model method is not applicable because 
accurate satellite parameters or orbit parameters are not available, or where 
the polynomial method is not applicable because uniformly distributed control 
points are difficult to obtain. 


Another advantage of the proposed method is that it can perform simultaneous 
rectification on a full-scale 200x1400 mm2 image without having to process 
the individual sub-images. This not only improves processing efficiency but 
also eliminates possible errors introduced by the joining process of recti- 
fied sub-images. 


It is generally not difficult to select more than 30 ground target points on 
a full-scale image. A sufficiently large number of uniformly distributed 
ground target points should be chosen at the two sides of a satellite 
panoramic image where large distortions occur. Simulation results show that 
this method can meet the accuracy re,i1irements for geometric processing 
satellite images. 


The authors wish to express thanks to Comrades Ju Tao, Guo Zhaozeng, and Li 
Lingzhi for their support in developing this method. 
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[Text] Abstract 


In this article, an improved method of extracting the poles of underground 
targets is presented. It is pointed out that in a noise environment, the 
eigenvalue method is more accurate in determining the number of poles; it is 
also shown that the eigenvalue method always yields smaller error of approxi- 
mation than classical methods. Experiments of extracting the poles of buried 
targets have been performed, from which several significant conclusions can 
be drawn. 


I. Introduction 


An effective approach for detecting underground targets is the transient 
electro-magnetic field method. This method has applications in many 
different areas such as hydrology, geology, resource exploration, detection 
of buried articles, mine detection, and anti-stealth technology. However, 

to use this method, it is necessary to solve the problem of target identifi- 
cation using the transient target response. Since the poles are invariant 
properties of a target, they can be used as effective signatures for target 
identification.1»2»3 In 1975, Rai Mittra et al. proposed the Prony method of 
extracting poles from transient pulse response. This method attracted a 
great deal of interest, but further study showed that in the presence of 
noise and measurement errors, the accuracy of this method is so poor that it 
may be of little practical use. Therefore, efforts are being made to explore 
new methods of pole extraction. 


II. Comparison of Square Approximation Error 
In 1978, Van Blaricum et al. proposed a new method of pole extraction’ where 


the problem is formulated as one of finding the smallest eigenvalue (Amin) 
and the corresponding eigenvector a of the matrix HTH: 


a™(H)"(H)a=4 (1) 
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where 


pf, | 
i i Ine, 
(H) > 
Tin Lime po 


I; is the sample value at time iAdt, and al = (ay, ‘**, G5); the vector a 
satisfies the condition: 


a 
> arZ’=0 (2) 


p=o 


where N is the number of poles. Once a; is determined, the poles in the 
Z-plane can be obtained from equation (2). 


Let € = (€55 Ej» °**s em) be the error vector, which can be expressed by the 
following relation: 


(H)a=e (3) 


If W(e) is the square approximation error, then we have: 


a’ (H)"(H)a= Soe? =W(e) (4) 


=0 


Different a vectors may yield different values of square approximation error; 
the a vector determined by the eigenvalue method is the eigenvector that 
corresponds to the smallest eigenvalue A,j, of the matrix HTH, i.e., 


a? (H)" (H)G@=Aein (5) 


Therefore, the problem chat needs to be investigated is the relative magni- 
tude between Api, and Wie). 


By arranging the eigenvalues of HTH in ascending order: 

O< hein ( =ho) SA SA yg SAy (6) 
then it is always possible to find a series of (N+l) orthogonal unit 
vectors Ej (i = 0, 1, °**, N) such that (u) = (E,, EE}, °***, Ey) is an 


orthogonal matrix, and satisfies the condition: 


\ 


Ami: (0) 
(u)™(H)* (HA) (u) | A, *, (7) 
Oo *, 


Consider a normalized vector a = (aps 4], ‘**', ay)! ; from matrix theory we 
can write: 
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a Dak or ad=(u)a (8) 


Therefore, 


li (ej=a' (A) (1)ea 


a'(uya’(H)'(H) (ua 


YA: 
-. (9) 


It can be seen from equation (9) that when all the eigenvalues are equal to 
Amin» we have the following inequality: 


I (€) DY) 4-0 A Sia (10) 


Since a is a normalized vector, it follows that 
Wie) >a (11) 


Thus, we have proved that the square approximation error of the eigenvalue 
method is always smaller than that of classical methods. 


III. Determining the Number of Poles 


The number of poles can be easily determined using the orthogonality 
property of the eigenvectors. 


1. In the noise-free case, let 

(p) = (/1)' (H) 
then the following procedure can be used to determine the number of poles: 
a. Construct a y x y matrix (¢) from the sample values. 


b. Find the eigenvalues of the matrix (¢) and determine the number of zero 
eigenvalues; if there are L zero eigenvalues, then the number of poles is 
V=y-L; if there is no zero eigenvalue, then the dimension of (¢) must be 
increased to form a new matrix (6). 


a 


c. When there is only one zero eigenvalue, i.e., when L = 1, then the 
corresponding eigenvector is the Prony vector a, and the poles are deter- 
mined along with V. 


2. In the case where noise is present, this method can still be applied. 
Let the noise sequence be zero mean with variance 02; since the noise 


NO 
Co 








sequence is uncorrelated in time and also uncorrelated with the signal 
sequence, the orthogonality property of the eigenvectors is preserved, 
Analysis has shown that the noise-free matrix (¢) and the noise-present 
matrix (¢') are related by the following expression: 


(d) =(¢) +¥0°l (12) 


where Il is the unit matrix. Therefore, the noise-free procedure can still be 
used, except that in this case, the eigenvalue sought is yo2, not zero. Once 
an eigenvalue of yo* is found, the corresponding eigenvector is the a vector. 


The above analysis shows that this method is very effective in treating the 
noise problem; once the two basic statistical parameters (mean and variance) 
of the noise sequence are known, the number of poles N and the a vector can be 
readily determined. 


IV. Experimental Procedure and Data Processing 
1. Experimental Set-Up and Procedure 


A block diagram of the experimental set-up is shown in Figure 1. The pulse 
width is t = 1.2 ns, the amplitude is 40 v, and the repetition frequency is 

10 KHz; the number of sample points is 256, the time duration is 20 ns, and 
the sample interval is 0.078 ns. The environment is a 100 cm x 100 cm x 60 cm 
sand box. The transient pulses are transmitted by a narrow-pulse generator; 
the transmitting and receiving antennas are quasi-TEM horn antennas. The two 
receiving antennas have identical characteristics; the received signal is 
amplified and stored in a micro-computer. The synchronization of the trans- 
mitted pulse and the sampling circuit is also controlled by the micro- 
computer. 
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Figure 1. Detection of Underground Targets Using Transient Pulses 
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2. Provessing of the Initial Waveform 


a. Elimination of the effect of ground reflection. In this experiment, a 
dual receiving antenna configuration is used so that the effect of ground 
reflection can be eliminated by differencing the two received signals. 


b. Noise processing. Noise reduction is accomplished by averaging several 
waveforms measured at the same antenna position. Any d.c. component is 
eliminated by adjusting the baseline. 


c. Gate switching. Gate switching technique is used to eliminate the 
direct coupling effects of the antenna and premature target response. 


d. Filtering of high-frequency components. Because of inconsistencies in 
sampling times between the two channels of the received signal, high- 
frequency "burrs" often appear. To solve this problem, a low-pass filter 
is used to process the waveform. 


The above pre-processing of the data significantly improves the signal to 
noise ratio, and consequently ensures the success of pole extraction. 


V. Conclusion 


Experiments of extracting the poles of various targets have been conducted; 
four of these targets are shown in Figure 2. The targets are buried in 
homogeneous dry sand at a depth of 3-10 cm; the receiving and transmitting 
antennas are separated by a distance of 40 cm, and the sampling interval is 
0.078 ns. The measured and recovered waveforms are shown in Figure 3, where 
the pole locations are also indicated. 
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Figure 2. Physical Properties of Four Different Targets 
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(- 0.00944, + 6.65178) 
t( ' ) (-~ 0.01199, + 0.85874) 






(- 0.01623, + 0.26459) 


_ 1 
Poles (X wT ) 


(~ 0.04658, => 1.10444 


a UAT) _o 0467] 1.42127 ) 
LV W\ _ oo 


y.% - 
“TV uso\/"  (- 0.01083, + 0.23539) 











(d) Target 4 


Figure 3. Poles of Buried Targets, Measured Waveforms and 
Recovered Waveforms 








On the basis of analysis and experimental results, the following conclusions 
can be drawn: 


1. The eigenvalue method always yields smaller error of approximation than 
classical methods. 


2. In the presence of noise, the eigenvalue method can still be used to 
determine the number of poles and to extract the poles as long as two basic 
parameters of the noise statistics are known. 


3. Changes in target depths (within 10 cm) has little effect on the poles; 
however, the measured results with respect to target depth are not presented 
here due to space limitation. 


4. The initial waveform must be pre-processed to facilitate pole extraction. 


5. The initial phases of conducting and non-conducting targets show clearly 
opposite sense (see Figure 3-a and Figure 3-b). 


6. It can be seen from Figure 3 that the larger the scattering object, the 
smaller the real part of the pole and the larger the imaginary part; con- 
versely, the smaller the scattering object, the larger the real part and the 
smaller the imaginary part. 
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a Three — Axis Flight Simulation Table 





40100067B Harbin HARBIN GONGYE DAXUE XUEBAO[JOURNAL OF HARBIN 
INSTITUTE OF TECHNOLOGY] in English No 3, June 90, pp 121-127 


Wang Mao Zhao Keding Liu Qinghe 
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Abstract In this paper , a pure digital control method for the outer axis of a three — axis 
flight simulation table is worked out and a state observer is used to obtain the 
acceleration signal .The method is simple and casy to be emplayed, and especially suitable 
for the electrohydraulic servo—system with fast response . The simulation results show 
that the system given in this paper has optimal three — order system performance . 


Key Words Electrohydraulic position servo system; state feedback; digital control . 


0 Introduction 


A flight simulation table is a kind of experimental equipment used for simulating the flying 
attitude of a tactical missile. Therefore, it is required to have the performance of high accuracy, 
fast response, less following time and little overshoot Now a days, the optimal linearization 
tested by an analog potentiometer can only reach + 0.02% (and the best position accuracy 
is + 4°32 ’). Besides, the deficiencies of parameters shifting and nonlinearity often occur in 
other analog components. Therefore, the analog control system today is difficult to satisfy the 
performance required. Furthermore, because the classical digital control is only restricted to sin- 
gle feedback loop system, digital design is very difficult, and especially there is no sufficient the- 
ory for a multi—loop feedback system, the hybrid dectrohydraulic servo system is widely 
used in almost all the control system design of the table. No doubt, this scheme is unfavorable 
to improving its accuracy and reliability. Of caurse, optimal control and adaptive control are 
convenient for digital implementation, but they will take a long time to be carried out in 
software and this can not satisfy the expected bandwidth of a system In addition, since the 
wordlength of a microprocessor is constant, a complex control algonthm will make large num- 
bers of accumulated errors. In order to avoid the above difficulties, a simple algorithm based 
on a state observer and state feedback is given in this paper. The digital simulation results indi- 
cate that its performance is satisfactory. 


1 Algorithm of the Controller 
In order to obtain acceleration signals and implement the system conveniently by 








@ microprocessor , a state observer is used to gain all the states in this design. The system 
block diagram is shown in Fig. 1. 
The state space representation of the system is given by 


X= AX + BU 


Y=CX (1) 
where { A, B, C } can be both controllable and observable. 
Assume the sampling period to be To, then the discretization of Eq. ( 1) is as follows: 


Mk+1) =FXAkK) + GU(k) 


¥{k+1) =CXk+1) (2) 


Te 


where F = ¢*% c-| e“"dt-B and Eq. (2) can be both controllable and observable. 
0 


Introducing state feedback, we obtain 


U(k) =r(k) —KXk) (3) 
Substituting Eq. (3) into Eq.(2), we have 


Mk+1) =FMk) + Glr(k) —KMk))=[F -GKJM k) +Gr(k) 


¥(k+1) =CXk+1) (4) 

From the linear system theory, we know that the poles of the system can be assigned to 

any expected position by selecting different feedback matrix K, and the introduction of state 

feedback will not destroy the controllablity of the system. Since Eq. (2) is observable, we can 

introduce a dosed loop state observer to obtain all the states . Generally , the observer has 
two forms: 


(1) Mk+1)==FRk) + GUk) 


Mk +1) = FXk+1) +L(Wk+1) -Wk+1)] 
¥(k+1) =CMk+1) | (5) 


where X[k) is the reconstructed state of X[k), and L is the gain matrix of the observer. 
(2) 


Mk+1) =FMk) +GU(k) +L[¥k) -Hk)] 


Hk +1) =CMk+1) (6) 
where L is the same as the above — mentioned. 
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Substituting X[k) into Eq.(3), we have 
U(k) =k) -KXk) (7) 
Substituting Eq. (7) into Eq.(2), we obtain 
Xk+1) =FXk) -GKXk) +Gr{k) 
Yk+1) =CXk+1) (8) 


In the continuous — time system, the separation theorem shows that the dynamics of the 
plant is independent of the dynamics of the observer, and vice versa. For a discrete — time sys- 
tem, there is also a same theorem. 

Let AX(k) =XUk) —ATk), then 


Xk +1) =( F-—GK) Xk) +GKAXk) +Grk) 

AX (k+1) =( F-LC) AX (k) (9) 
Rewrite Eq. (9) in the vector form as follows 
T Mk+1)) F-GK GK Xk) G) 

e + nk) (10) 
| AM k +1) 0 F-LC} L AMk) 0. 


— 











From the matrix theory, the characteristic polynomail of the discrete transfer matrix of 
Eq. (10) is given by Eq. (11) 


det (A 1— F +GK) (Al- F+LC)] = det{Al—F+GK] da (AI-F+LC]=0 (11) 


So the poles of the plant and the observer can be assigned to the expected position by 
separately selecting K and L. Consequently the controller and the observer can be designed 
separately. 

The stability of the system is discussed as follows. From Eq.(2), so long as the whole 
cigenvalues of the transfer matrix are inside the unit circle , Eq. ( 10) is stable. When the system 
comes into a steady-state, Eq. ( 12) is obtained. 


k-@ 

Eq. (12) shows when k is sufficiently large, the state error of the observer and the plant is 
zero, ie, the states estimated by the observer are the real states of the plant 

Generally, the whole state initial values of the plant are impossible to be known, so, 
when X(0) # 0, at the beginning of the system transient procedure the observed states are 
not the actual states of the plant; and when X[0) =0, the observed states are the actual ones 
of the plant at the beginning. Analysing Eq.(6), we find that if X{0) =O, then the observer 
matrix L is useless and an unclosed —loop observer will satisfy the requirements. 


lim 4X{k) =lim(F-LC)}* 4X0) =0 (12) 
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2 Pure Digital Control for the Outer Axis of a Three—Axis Flight 
Simulation Table 


The electrohydraulic servo system used in the outer axis of a three — axis table can be ex- 
pressed approximately as follows. 
G(s) = s(s*+2éa,s +m ) 





where 


D,, = 800ml/rad — displacement of the motor 
w» = 207.35rad/s — characteristic angular frequency of the motor 
& =(.1 — damping ratio of the motor 
' K— gain of amplifier (determined by speed errors) 
The following problems will be considered in design. 
2.1 Selecting Expected Poles 
From Eq.( 13), we know that the system is of three—order. Then, we can make the ex- 
pected poles of the plant equal to the dosed—loop poles of a three — order optimal system. 
The expected system is expressed as follows. 


3 


or 
Cds) = S35 70,s9+ las tos 





(14) 


where w ,= K,, and K, is the open loop gain of the system. 
Here, we take K,= SO and substitute it into Eq.( 14}, and then we obtain the expected 
poles of the plant as follows. 


A,= — 106.0639, 2,,;= —47+ j85 


In order to obtain the actual states of the plant as fast as possible, the expected poles of 
the observer should be as far away from the Im axis in [s] plane as possible, but in practice 
it is limited by the oscillating frequency of the dock Generally its real parts are taken as 
(2 ~ 3) times as the real parts of the expected poles of the plant Here, we take 
A,=2A,= —212, A,,;=2Re{A,,]+ (0 


2.2 Selecting the Sampling Period 


The sampling period T, must satisfy the sampling theorem. When 7, is too large, some- 
times we can obtain an expected robust discretization system, but after state feedback being in- 
troduced the dosed-loop system will be an unstable positive feedback system ( the 
unclosed —loop system must be stable ). When T, is too little, the accumulated error will 
make the system unstable and reduce its accuracy. 

So, from the viewpoint of engineering, T, will be selected as follows. 
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w,> max(5~10) x(a, or the bandwidth of the system) 
In this paper, the characteristic frequency of the observer is larger than that of the ex- 
pected system, and it has decisive influence on the system So w, should be the characteristic 
frequency of the observer. Then 





Here, T, is selected as 0.0005s. 


2.3 Construction of the Control System 


The implementation block diagram is shown in Fig 2 The controller and the observer in- 
side the dosed dotted line are implemented by a microprocessor. Because the introduction of 
state feedback damps the output of the system, we introduce a compensator into the input 
pass. Its parameters are determined as follows: 

From the viewpoint of the continuous system and assuming that there is no observer, the 
system shown in Fig 2 is equivalent to the system shown in Fig 3. Because k, is the feedback 
gain of state x, and considering that Y=c , x,,the output feedback gain in Fig 3 should be 
k,/c,. For the unit step input, the output of the system is 








l c l 
Ys) = = = 
(sG(s)+—+-c,) 5 sG{s)+k,'s 
i 
Then y(co) hm s. Hs) =¢,/k, 


So, in order to ensure y(oo ) = 1, the block ? in Fig 3 should be replaced by an amplifer 
with a gain k,/c ,. 


3 Simulation Results 


Based on Eq.( 13), digital simulation has been done. The step and sine response of the 
plant with zero and little and large initial states are discussed respectively. When the state ini- 
tial values of the plant are zero, the system is equal to one without an observer. Its response is 
shown in Fig 4. The unit step response of the plant with initial states is shown in Fig 5. 
When the state initial values of the plant is not zero, a non-normal balance point occurs in 
the position following time .Comparing Fig. 4 with Fig. 5, we find that the position following 
time is about 0.05s and 0.1s respectively for the system without and with an observer. The sin 
(79t) response (for the gain frequency bandwidth test) and the sin (69t) response ( for the phase 
frequency banduidth test ) of the system are shown in Fig . 6a and Fig . 6b respectively. 


4 Conclusions 


The pure digital control system for the outer axis of a three — axis flight simulation table is 








designed in this paper. It is simple and convenient for implementation. Simulation results show 
that although a non-normal balance point occurs in the system with an observer, the scheme 
given in this paper can ensure perfect dynamic performance for the three—axis table and suit 
for a most dectrohydraulic position servo system. 












Fk +1) =FX(k) +Gu(k) | 


> +L (¥(k) —¥(k)) 
¥(k+1) =(X(k +1) 
Ui kK) =rik) —§ 





Fig . 1 A discrete-time control system including an observer 
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Fig .4 The unit step response of the plant Fig. 5 The unit step response cf the plant 
with a zero initial value with a little initial value 
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Fig . 6 The sinewave response of the plant with an initial value 
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COMPUTERS 


A Hybrid Learning Method for Multilayer 
Perceptrons 


une 


40100067A Harbin HARBTY GONGYE D UF SHEBAOT JOURNAL OF HARBIN INSTITUTE OF 


OF TECHNOLOGY] in Zhao Meide Huang Wenku ‘Hong Jiarong 
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Abstract A Newton learning method for a neural network of multilayer perceptrons is pro- 
posed in this paper. Furthermore, a hybrid learning method is legitimately developed in 
combination of the backpropagation method proposed by Rurrelhart et al with the 
Newton learning method Finally, the hybrid learning algorithm is compared with the 
backpropagation algorithm by some illustrations, and the results show that this hybrid 
learning algorithm has the characteristics of rapid convergence. 

Key Words: Computer; Newton learning method; neural network; multilayer perceptron 


0 Introduction 


Multilayer perceptrons are feed — forward nets with one or more layers of nodes between 
the input and output nodes. These additional layers contain hidden units or nodes that are 
not directly connected to both the input and output nodes. Multilayer perceptrons overcome 
many of the limitations of single layer perceptrons ™ , but were generally not used in the 
past because effective training algorithms were not available. This has changed with the devel- 
opment of new training algorithm, i ¢, the backpropagation algorithm ™. Since its 
emergence, they have been quite successfully used for a wide vanety of tasks such as informa- 
tion processing, coding and pattern recognition °“? . Unfortunately, in many applications, 
the number of iteration required before convergence can be large and there exists local minima 
problem. Some modifications to the backpropagation algorithm were given by Vogl et al™ in 
order to accelerate the convergence of the backpropagation algorithm. But we found that 
these modifications cannot accelerate the convergence by many numerical simulations. In Sec- 
tion 2 of the paper the backpropagation algorithm and some modifications to it are briefly des- 
‘ ctibed. In Section 3 the recursive formulas for calculating Hessian matrix of multilayer 
perceptrons are derived and a Newton leaming algonthm is given. A hybrid learning algorithm 
for multilayer perceptrons is proposed in Section 4, and:in Section 5 two examples are 
given to compare the two methods. 


1 The Backpropagation Method and Its Modification 
1.1 Model of Neural Nets 


4O 





In general. a neural! network consists of N nodes ( processors), cach of which is con- 
nected to all the other nodes. In the present work, however, only layered feedforward networks 
are used .Use of such a layered structure sacrifices a certain amount of parallelism since only 
processors in the same layer operate in parallel but adjacent layers compute scraily. Although 
some parallelism is lost, these networks are popular because there exists a sophisticated learn- 
ing scheme for training the metwork. A typical network using “¢ layered structure is shown in. 
Fig. 1. It consists of an input layer,an output layer, and hidden layers in between. The nodes 
in the hidden layer are necessary to implement nonlinear mappings between the input and 
Output patterns. | 

Associated with each node is a threshold value. The output of each node is_ determined 
by the level of the input signal in relation to the threshold value. Signals can propagate 
through the network only between adjacent layers and only in a single direction. These sig- 
nals are modified by the connection strengths (.weights) between the nodes. Let 

6, , threshold of the ith node in the kth layer, 

W, ,.. weight betweenj th node in the (k—1)th layer to ith node in kth layer, 

net, , input to the ith node in the kth layer, 

O,, Output of the ith node in the kth layer. 

The input to a node is given by 


nel, =) [ Win, 2-1] + 0, (1) 
j 
The summation index j extends over all nodes in the (k —1)th layer. The output of a node is 
given as 
O; a-1 =j (net, 4) (2) 
where f is a sigmoid function given by 


j (net x) = a (3) 





1.2 The Backpropagation Method and Its Modification 


In order to describe conveniently, here we use vector notation. The basic iteration formu- 
la of the backpropagation algorithm is : 


W**! = W*-n, VE(W*) +0, (WY Wr") (4) 


where W is a column vector which involves all the weights and the thresholds, W* represents 
the column vector of W in kth iteration, E( W) is the total error, ie , the sum of the squares 
‘f the error of the output signal , V E( W) is the gradient coli nn vector of H W) with respect 
to W, n, is a learning constant , and 2, is a momentum factor. 

In the backpropagation algorithm, n, and a, are preset constants. Lea E(W) =Ej W), 
where E,( W) is the error function under the pth input pattems, Le, 
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E,(W) => E (4y- Op.) (4) 


where the summation index j extends over all nodes in the output layer. _, is the desired out- 
put ( provided by a teacher via the input-output samples) at the ith output node and O, ,, is 
the actual output at the ith output node Ths element of V E( W) is: 


ai. = 9 1 0p (5) 

5 ix is the error signal for unit i of the kth layer under the pth input samples. 
; j (net yin (GeO 4)  jork=n (6a) 
me j net yi.) 2 (5 pyre W, +1) otherwise ( 6b) 


where the index s indicates summation over al! node: °° the kth laver. The detailed description 
of the backpropagation algorithm is given in [2] anc ~; 

Some modifications to the backpropagauor. aigoritnm was given in [5]. These modifica- 
tions are as follows: Instead of updating the network weights after each pattern is presented 
tc the network, the network is updated only after the entire repertoire of patterns to be learn- 
ed has been presented to the network; Instead of keeping » constant, it is varied 
dynamically so that the algorithm utilizes a near—optimum » _; the momentum factor a is 
set to zero when, an signified by a failure of a step to reduce the cotal error, the information 
inherent in prior steps is more likely to be misleading than beneficial. Only after the network 
takes a useful step, ie, one that reduces the iotal error, does a again assume a non-zero 
value. 

Here, we propose a new scheme for selecting x, and n,, they are selected such that 


O%, | = min E(W*-n V E(W*) +o W*-W’)) 
a,” 
Theoretically we have solved the problem of how to select « and » in the backpropagation 
algorithm. However, by many numerical simulations we find that both the modifications to 
the backpropagation algorithm given in~[5] and above scheme cannot accelerate the 


convergence of the backpropagation method. In the next section we propose a new learning 
method for multilayer perceptrons. 


2 The Newton Learning Method for Multilayer Perceptrons 


In order to solve above problems, we propose a Newton learning method for multilayer 
perceptrons. 


2.1 The Hessian Matrix of Multilayer Perceptrons 


The Hessian matrix of the network under the pth input patterns is denoted as H, Hence 
we have 


4o 











H,=9°E, 
where V*E, is a second order derivative matrix of E, with respect to W. For convenience, we in- 


troduce the following notation: 














H’ oF = @E, 
tut OW. OW, |. (7a) - 
05 
Yud,e = Gnet. ( 7b) 
, a5 , 
iJ.ee@ Onet ») ( 7c) 
: _ 00,,., 
Wr Get (7d) 
Since H, is a symmetric matrix, we can discuss it in two cases : 
1) T he diagonal elements of H, a; A 7 
From (5), we have, 
0 
Ba, 7 HO 
By (7a), we obtain 
- VE, 66 
Aoi = OW, , OW, , = On, + -1 WW, 
since 
By, _ 00, Onet, ;., 
ow, Onet »;,, OW... 
According to(1) and (7b), we can obtain 
OW, ., plirnlfi.j.r 
Hence 
Ay) = — “phir U5) p=1 Vivy,r 
( 8a) 


2) The non —diagonal elements of H, Ay (r >q) 
From (7a) and (5), we have 


O ( -4,, ,On.,-;) nets, 





Aus ~ Onet, ;4 , OW, 1.4 
Onet , , ‘ 06, , On, r-1 
Zs O r~ + 06, ’ Onet 
OW, 1.4 ( = Onet »,,4 ™ en 
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Acoording to(1), (7c)and(7d), we can obtain 


HS, alin Os 1, ¢-1 (O..-) Cure +35, Day t.e) ( 8b) 


where the recursive formulas for calculating y, ,,,{,,,¢ and 9), ame given in Appendix. 


2.2 The Newton Learning Method for Maltilayer Perceptrons 


According to the Hessian matrix derived above, we can obtain the basic iteration formula 
of the Newton learning method for multilayer perceptrons as follows: 


w' = W-y,H'G, (9) 
where H,=) H,=) VE, (W), and 
r : 
G.-T VE(W) 


The index p represents summation over all input-output samples used for training. y, is a 
constant which is called a Newton lea ning constant in the paper. 

The Newton learning method has been implemented on a Sun3 Workstation in Clanguage. 
We find by many illustrations that this algorithm cannot converge to a solution. There are pos- 
sibly several reasons for it through our analysis: 1) Although the Newton algorithm can 
quadratically converge near the exact solutions, it is very difficult to find : better initial guess; 
2)After several iterations, the matrix H will become singular ;3)We cannot ensure that for each 
iteration the quadratic model can adequately approximate the objective function E. 

To apply the Newton leaming method proposed in the paper, we propose a hybrid learn- 
ing method which combines the backpropagation method with the Newton learning method. 
The hybrid learning method has the advantages of the two methods. The illustrations given in 
Section 5 show that this new learning algorithm has a surprising speed of convergence. 


3 The Hybrid Learning Method for Multilayer Perceptrons 


Since the backpropagation method is a steepest decent method, it can globally converge to 
a solution. The Newton leaming method can quadratically converge in a local region of 
the solution. Therefore, the hybrid learning method can rapidly globally converge. The 
basic idea of the hybrid leaming method is: after each input pattern is presented to a net- 
_ work, the weights and thresholds are updated-by the backpropagation algonthm; after entire 
input patterns are presented to the network, the weights and thresholds are updated by the 
Newton learning method. The detailed description of the algorithm is as follows. 
Step! Initialize weights and thresholds (set all weights and node thresholds to small ran- 
dom values ). 
Step2 According to (1) and (2), calculate the actual outputs of cacn node under in. 
pth input pattern in forward direction. 


4h 











Sicp3 Calcuiate che error signals of the ouipui noues in cerms of (6a) and update the 
weights and thresholds of the output layer according to ( 3) 

Step4 Calculate the , ,, of the output layer in terms of (A 1a) and calculate the corre- 
sponding elements of the matrix H, by (8a). 

Step5 Calculate (, ,,, andg,,, ,in terms of (A.3) and (A 4), then calculate the corre- 
sponding clements of the matnx H, by ( 8b). 


»), Of the penultimate nodes, and calculate the corresponding 
cements of the matnx H, by ( 8a). 
Step7 Calculate the error signals of the penultimate nodes in terms of (6b), and update 
the weights and thresholds of the’ penultimate layer according to (3). 
Step8 Repeat Siep 5, Step 6 and Step 7, until the entire elements of matnx H, is calcu- 
lated. 
Step9 Repeat Step 2—Step 8, until all input patterns are presented to the network. 
.Stepl10 Decide whether or not matrix H has an inverse matnx If it has, then go to Step 
11, otherwise go to Step 12. 
Step|1 Update all weights and thresholds of the network according to (9). 
Step12 If the maximum of the differences between the actual output and the desired out- 
put of the output layer for all input patterns is less than preset value, then termi- 
nate, otherwise go to Step 2 


Step6 By(A1 )calculate the » 


4 Illustrations and Comparison 


The backpropagation algonthm and the hybnd algonthm have been implemented on 
Sun3 Workstation in C and FORTRAN languages. Two examples are given as follows. 


4.1 The XOR Problem 


The XOR problem is a difficult classical problem, which can’t be solved by single layer 
perceptrons, Since it is a simplest nonlinear seperatable problem. The XOR problem, ie, Ex- 
dusive OR problem is the problem that a neural net can generate the following outputs for the 
following input samples by training. 


Input Pattem Output 


0 0 0 
l 1 0 
0 1 l 
l 0 l 


Fig. 2 and Fig. 3 show the solution solved by the backpropagation algonthm and the hybrid 
algorithm respectively. The solution in Fig. 2 is reached after 532 sweeps through the four 
stimulus patterns with learning rate of » =0.S.and momentum factor of 2 =0.9. The solution 
in Fig. 3 is reached after 39 sweeps with » =0.5, 2 =0.9 and y =0.11. 

Through many numerical simulations we find the parameter » and a effects on the num 
ber of iterations before convergence for the two methods. Fig .4 and Fig. 5 show these effects 
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respectively. The Newton learning constant y affects the number of itcrations before 
convergence for the hybrid algorithm Fig .6 shows this effect y is now varied between 0.005 
and 0.02 by experience. 

In numerical simulation we also find that the Newton updation for the waghts and 
thresholds only ocou in the first several sweeps in the hybrid algorithm Table | is ob- 
tained by assuming that afler 10 sweeps the Newton updating stops. 


Table L Set a=0.95 and y=0.011 





























: Backpropagation | Hybnd Leaming 
sweep | CPU tre | sweep /|CPU tre 
0.1 | 2736 | 246 127 6 
'o2 | 1011 87 124 6 
cs SBI 51 6 | 4 
Qe 383 35 29 2 
| 0.5 | 274 25 28 2 




















BP — Backpropagation HL —Hyurx Learning 


4.2 The Parity Problem 


The parity problem is a problem which can computer the parity of any input pattem, 1 
e, the output required is | if the input pattern contains an odd number of is and 0 other- 
wise. This is a very difficult problem because the most similar patterns require different an- 
swers. The XOR problem is a parity problem with input patterns of size two. Here we have 
done a parity problem with patterns of size 3. Fig. 7 and Fig. 8 show the actua! solutions with 
three input lines and three hidden units obtained by the backpropagation algorithm and the hy- 
brid algorithm respectively. Let » =0.5, « =0.9, and y =0.006. The solutions in Fig 7 and 
Fig. 8 are reached after 1441 and 52 sweeps through the eight stimulus patterns respectively. 


5 Conclusions 


In this paper, a Newton learning algonthm for multilayer perceptrons is derived and a hy- 
brid learning algorithm is develop.d in combination of the backpropagation algonthm with 
thc Newton learning algorithm . The results of numerical simulations given in the pape. siinw 
that the hybrid algoxthm has much faster convergence than the backpropagation algorithm 
The future study of the project is to generalize this new aigorunm: to high oraet ncurses nd 
apply it to solving practical problems. 
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Appendix 
The Recursive Calculation For y, ,, Ciicg ANd iin 


1. The Recursive Calculation For ;,, 
We derive the formula of y; ;, in two cases as follows. 


(i) Layer ris an output layer 
From (6a), we have 


911 =(ds Op, ’ (net ,,., ) 


By (7b), we can obtain. 


» = {° jor i#j 
UG A net ys Ms Oi.) LG’ (net ys.) F, otherwise —_ 
(ii) Layer r is a hidden layer 
From (6b), we have 


551.7 = j (net,,,)> 8) ret Ww, el,r+l 
J 


(1) For i#j, 
Mtr =f ne VL Mires * Gnet y, , 
By (7c), we obtain , 
Yje =f (net s.r) d Watirst Cay rete) ( Alb) 


(2) Fori=:, 


Mj, =j “(net 51.) ¥( Oy srt Wieisrot) +j (net yi) L(W,. 1 Seutsretied 


2. The Recursive Cuiculation For, ,. ;.q ( Alc) 


According to (7c), we have. 





06 06 Onet 
= pir = pi? 21a? 
Cisne Onet L( Onet 


pr.qri 


5 60 


Pj.@ 
= —-W r+ net 
inet y).4 ~~ O04 Omety gay TUM) (a2) 








Basa {Ber + br rma 


bs.0.0,¢08 , jor r>q+l 
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¢ j (net ys) (W, ptn0d Ve.a.08 2 jor r= q+l1 
i,j.t aa 
j'(nét »,,, LW, eJert Ci arger) » Jor r>qgtl (A. 3) 


3. The Recursive Calculation For 9; ; ,., 


(i) For r=q, 
By (7d), we obtain , 


60,,,  _ ee 





Pi, »?,? as . ° 
Onet yy, (J {net,,,) jor i=j (A, 4a) 
(ui)Forr >q 
According to (7d) and by chain rule , we have 
_60,, 00,5. r Onet ys .¢>1 
9, Loe . Te > ae Onet p ,. 9 énet _ 


By (A. 2), we obtain 


P, j,r, ¢ = j {net 5) ). d ( W, 8.J, e+! P,, s,? Pay (A4b) 


c2 ce oe aa 


am >) > &® C>| 


— 
(k=1)t yer) CBD OD 








Se OG 


Fig. 1 The structure of a feed-forward net 
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